Time delay in the Kuramoto model with bimodal frequency distribution 
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We investigate the effects of a time-delayed all-to-all coupling scheme in a large population of 
oscillators with natural frequencies following a bimodal distribution. The regions of parameter 
space corresponding to synchronized and incoherent solutions are obtained both numerically and 
analytically for particular frequency distributions. In particular we find that bimodality introduces 
a new time scale that results in a quasiperiodic disposition of the regions of incoherence. 



The Kuramoto model is presumably the most suc- 
cessful attempt to study macroscopic synchronization 
phenomena arising in large heteroggneous populations of 
interacting self-oscillatory units PH^il^l. Kuramoto, mo- 
tivated by Winfree's work on biological oscillators 0, 
showed that the dynamics of an ensemble of N weakly 
interacting limit-cycle oscillators can be treated consid- 
ering simply the oscillator phases [9i, . . . ,9iq). In this 
paper we study the Kuramoto model with delayed inter- 
actions M 



K ^ 
N ^ 



sin[0,(t)-0,(t-r)]+e»(i); i = 



(1) 

where heterogeneity is established considering a certain 
distribution giio) of the natural frequencies w^. The 
terms f i (i) represent uncorrelated zero- mean white noise 
processes, = 0, (6(00 (^')> = Wb,,b{t - t'). 

In absence of time delay (r = 0) and for large iV, as 
the coupling strength K exceeds a critical threshold 
the model Q shows drastically different transitions to 
collective synchronization, depending on the shape of 
9{^) [lIlllllllH- For a strictly unimodal distribution 
the transition occurs between a totally incoherent state 
and a partially synchronized state. In contrast, symmet- 
ric bimodal distributions give rise to hysteresis and/or 
a transition to a time dependent state composed of two 
clusters 

The interactions in ensembles of coupled oscillators 
have been traditionally considered to be instantaneous, 
an assumption that considerably simplifies the analysis 
of such systems. However, the study of phase oscilla- 
tors with time-delayed coupling j23| is receiving interest 
since a number of theoretical studies show that time delay 
may considerably affect the synchronization phenomena, 
typically leading^o multistability of many synchronous 
states (see e.g. |3 and references therein). In par- 



ticular, the Kuramoto model with unimodal frequency 
distribution has been generalized to allow time-delayed 
interactions in |^ ^| . Additionally, phase models with 
time delay have successfully explained synchronization 
between plasmodial oscillators jXj and in semiconduc- 
tor laser arrays Recent studies also demonstrate 
that time delay may be a useful synchronization-control 
mechanism in large oscillatory populations 13]. 

In this paper we investigate the effects of a bimodal 
frequency distribution on the Kuramoto model JQ) with 
time delay r. For r = and assuming that giio) is sym- 
metric (centered at w with twin peaks of width 7 at both 
sides), one may always transform to a rotating frame, 
such that the model (Q) is symmetric under the reflection: 
(Bi,u)i) [~Qi,~bJi). Time delay generally breaks such 
symmetry, except for the specific values t = Tn = wk I 
(n e N). Thus, at r = r„ we expect to recover the typi- 
cal reflection-symmetric stationary wave solutions found 
in the Kuramoto model with bimodal frequency distri- 
bution but without delay 0- For general r values the 
breaking of the reflection symmetry should give rise to 
the structures already observed for models without time 
delay but with either asymmetric bimodal frequency dis- 
tributions [1^ or asymmetric coupling functions |l5j| . 

We begin our analysis by considering the noise-free 
case, D = 0, and a frequency distribution that consists 
of two infinitely sharp peaks (i.e. 7 = 0): 



g{uj) = {8{lo - uji) + S{u! - uJii)]/2 



(2) 



with uji < ujii. Even for this simple choice the results 
are illustrative and far from trivial. 

Fully synchronized states.- Let us first investigate the 
existence and stability of synchronized solutions of 
consisting of two clusters of identical oscillators that ro- 
tate uniformly with angular velocity and phase differ- 
ence P 



(3) 
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FIG. 1: (color online) Frequency Q, phase difference P and 
leading Lyapunov exponents A (red), fj. (blue), for the stable 
synchronized solutions @ as a function of K for Auj = n/15, 
ijj — 77/2 (dotted line, upper panels). 



dental equations for f2 and /? 

n = Co - Ksin{nT)[l + cosf3]/2, (4) 
sin/3 = Atj/[if cos(rJr)], (5) 




FIG. 2: (color online) Stability boundaries for Synchroniza- 
tion and Incoherence. Theoretical boundaries for (full) syn- 
chronization are obtained from Egs. 16171 . Disintegration 
(A — > 0~) and unlocking —> 0~) are indicated with red- 
dashed and blue-solid lines, respectively. Black-thin curves; 
Boundaries of stable incoherence [Eq. Il2ll 1. Shaded areas: 
Incoherence regions obtained by numerical integration" of 
Eqs. Q. Dark-shaded regions: bistability between incoher- 
ence and synchronization. 

"Using a fifth order Adams-Bashforth-Moulton scheme with dt = 
t/20, except for r < 1, where dt = 0.05 , Af = 24, <Ii = tt/2. 



where Auj = ojjj — uji and (D = (w/ -|- uJii)/2. These 
equations have multiple solutions for given uj, Aw, 
and r [see Fig. Q]. 

The linear stability analysis of the solutions (jSJ yields 
two (y — l)-fold degenerate eigenvalues 

A± = -K[cos{VIt) + cos{Q.T ± /3)]/2, (6) 

and a set of eigenvalues {^.i} which are determined by 
the transcendental equation 

(7) 

where c± = cos(r2r ± /3), and c = cos(r2T). The eigenval- 
ues A+ , A_ govern the stability within the two individual 
clusters, each with oscillators. The eigenvalues (dis- 
carding the trivial solution ^0 = 0) are related to the sta- 
bility of the frequency locking between the two clusters. 
We denote the leading Lyapunov exponents in each sub- 
set by A = max{A+, A_} and n = max{Re(/.ti)} (i 7^ 0). 

Figure n displays numerical solutions of Eqs. I@J, lO, 
© and 10 for specific values of the parameters t, Acj 
and uj. The upper panels show that for r = r„ = wk /uj 
symmetry-related solutions with frequencies ^{K, t„) = 
uj±C,{K, Tn) arise in pairs around the "central" frequency 
u) (which is the synchronization frequency for r = 0) 
[Fig.n(a,c)]. For positive K, the central solution is sta- 
ble (unstable) for even (odd) values of n, and vice versa 
for negative K. For intermediate values, r„ < r < t^+i^ 
the stable solutions continuously vary between these two 
patterns, all of them approaching 17 = as r increases. 



This effect, common to time-delayed interacting oscilla- 
tory systems, is known as frequency suppression |16|. 

The growing number of synchronized solutions as K 
and r are increased (see Fig. was already reported by 
Schuster and Wagner for a system of two coupled os- 
cillators, corresponding to model 10 for = 2 without 
the self-coupling term. They also have shown that as K is 
increased, stable solutions alternately appear with small 
(/3 « 0) and large (/3 « tt) phase differences- this is also 
the case in the presence of self-coupling |23|. However, 
for a population — due to the existence of the destabiliz- 
ing modes linked to © — all the stable solutions are of 
the in-phase type, i.e., < /3 < 7r/2 [see Fig. cen- 
tral panels]. These solutions appear at saddle- node bi- 
furcations where either A or vanishes [see Fig. lower 
panels]. Increasing K, the Lyapunov exponent A is sta- 
bilized (A cx —K) whereas /i first decreases steeply and 
then abruptly changes its tendency (at the point where 
a pair of complex conjugate eigenvalues become the 
leading ones), and asymptotically approaches zero. 

Figure [3 shows the boundaries of the synchronous 
states Q in the t — K plane. For a unimodal distri- 
bution (Acj = 0) the regions without stable synchronous 
solutions are disconnected from each other, and the in- 
stability is always via the disintegration of the cluster 
[red-dashed lines, Fig.|2Ia)]. By contrast, for a bimodal 
distribution (Aw > 0) the stability boundaries detach 
from the K — Q axis and we obtain two separated con- 
tinuous curves of marginal stability for K < Q and K > Q. 
Note that as Aw is increased the instability of the phase 
locked state occurs mostly via unlocking of the two clus- 



3 



ters. 

Incoherent state.- Next we address an important type 
of solution of Eqs. QJ, so-called incoherence, in which 
the system is in a completely phase-disordered state. To 
study the stability of the incoherent state we introduce 
the complex order parameter Re'''^ = N~^Yl!j=i^^^^i 
that measures the degree of "phase coherence in the 
system. This permits to write the system Q in terms of 
the time-delayed mean field quantities R and ij} 



e,{t) =L0,- KR{t - r) sm[e^{t) ~ i>{t - r)] + i,{t). (8) 

Considering the limit N ^ oo, we drop the indices 
and introduce the probability density for the phases 
p{6,t,u;) ^1§\. Then p obeys the Fokker-Planck equation 



dp/dt = -d{pv)/d9 + Dd^p/d9^ (where v 



-KR{t^ 



T)sin[6'(<) — — r)]), for which the incoherent state 
Pq = (27r)^^ is always a trivial stationary solution. Lin- 
earizing the Fokker-Planck equation about po one finds 
that the stability of incoherence is determined by the 
eigenvalues A satisfying 



e — 



A + D 



-dhJ — 1. 



lUJ 



(9) 



This equation with the distribution Q and D — Q gives 
(after setting A = — iJ7c, to obtain the instability thresh- 
old) the collection of critical curves 

if« = (-2)'(a;,-17W)(^,,-f7W)/(cS-f](0),(io) 
17« = {\/2 + l)n/T^ (11) 

[l integer), and two trivial ones Kc = (for ~ toi and 
ujii). This leads us to the bounds for stable incoherence 

if+(r) = min{i^(')(r)|if«(r)>0}, (12a) 
K-{t) = max{i^(')(T)|ii:^')(T) < 0}. (12b) 

Figure 12 (a) shows the regions of incoherence obtained 
from Eq. (|12|l for a unimodal distribution (Aw = 0) 0. 
They consist of a periodic disposition of disconnected 
tent-shaped regions centered at t = t„ , of height Kc = 
(— l)""''^7r/T„ and width ti (at K = 0), as can be seen by 
setting u! = uji = uju in Eq. (^UJ. 

In contrast, the regions of incoherence for the bimodal 
frequency distribution Q are organized quasiperiodi- 
cally, with a decaying envelope of period 27r/Au; [see 
Fig-E^b-d)]. This can be understood noticing that, for 
\K\ w 0, such regions correspond to the overlapping 
of the incoherence regions of two independent unimodal 
populations (with uj = ui and loh). Thus, the new dispo- 
sition of the regions of incoherence is quasiperiodic with 
frequencies loi and lj//, and it can be determined through 
a function 

/(t) = cos(w/r) cos(ci;/7T) cx cos(2[I'r)-f cos(Aa;T), (13) 

which takes positive values exactly where incoherence is 
stable for some K ^ 0. 
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FIG. 3: (color online) /S.U — K diagram showing the regions of 
stable incoherence for a bi-Lorentzian frequency distribution 
of width 7 = 0.1 centered ui = tx 12. Solid-black (7 = 0.1) 
and brown-dashed (7 = 0) curves: analytical boundaries ob- 
tained from Eqs. I|14^ and l|12|l . respectively. Shaded areas: 
Regions of incoherence obtained by numerical integration of 
Eqs. CJ, 7 = 0.1, D = 0, = 2000. (a) Green curves 
(r = 0.1): symmetry breaking which results in two Hopf bi- 
furcations (the solid green curve becomes the new boundary). 
To the right of the red-dashed line at Aoj — 2'y/\/i [only 
shown in panel (a)] the bi-Lorentzian distribution is bimodal. 



Effect of frequency diversity (and noise).- For a uni- 
modal frequency distribution the inclusion of diversity 
(or noise) does not alter significantly the scenario already 
captured for identical oscillators roj , but for the bimodal 
distribution it has a more intricate effect. 

We restrict our analysis to the bi-Lorentzian distri- 
bution, but a qualitatively similar scenario is expected 
for other bimodal distributions. Specifically, we take 
giu:) = 7/(2^)([(c. - ioif + 72]-! + [(c. - uou f + 7']-^) 
which, for Alj > 27/\/3, becomes bimodal. Then the 
linear stability of the incoherent state, determined via 
Eq. lO, yields the transcendental equation 



A+ = 



(14) 

D and 7 appear at linear order always in the combina- 
tion (Z? -f 7), hence both effects can be studied together 
restricting to the case D = 0. 

With the help of Fig. JSJ^a), let us first recall the 
main results for t = [9|. For small separation of 
the peaks (27/\/3 < Aw < 27), incoherence becomes 
unstable subcritically in a steady-state bifurcation at 
= 2[(Aw/2)2 + 72]/7 |2||. In this case the system 
exhibits a single cluster of oscillators locked to the fre- 
quency r2c = ^- If Aw > 27, a degenerate Hopf bifur- 
cation at — 47 precedes the previous instability. At 
two symmetric clusters at both sides of the central 
frequency w appear simultaneously. The loci of both bi- 
furcations coalesce at T (a double-zero eigenvalue point 
in the co-rotating frame) . 

A similar scenario is recovered periodically at the val- 
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FIG. 4: (color online) t — K diagram with the stable regions 
of the incoherent state for the bi-Lorentzian frequency distri- 
bution for 7 = [dashed brown, from Eq. Il4ll 1. and [from 
Eq. III2I1 I for 7 — 0.05 (solid green) and 7 = 0.1 (solid black). 
Shaded area: Incoherence regions obtained by numerical in- 
tegration of Eq. 13, (7 = 0.1, D = 0, iV = 2000). 

ues T = T„, i.e., when system has reflection symmetry. 
The steady-state bifurcation occurs at the same values of 
K ^ but with a change of sign for odd n 

K!{r = Tn) - (-1)" 2[(Ac^/2)2 + 7^1/7. (15) 

Beyond the codimension-two point T, located on the 
curve iff, at Aw = 27-^/ (1 -t- 7r„)/(l — 7r„), this insta- 
bility is preceded by a degenerate Hopf bifurcation [see 
Fig. El^bjC)]. The latter expression diverges at r = 7"^ 
and hence, above a certain value of n (n > 4 for the pa- 
rameters in Fig. 131), another degenerate Hopf bifurcation, 
denoted by in Fig. E^d) (dash-dotted line), crosses 
the line and prevents incoherence to be stable for 
larger K H. 

A general value of the time delay, r ^ t„, changes 
this scenario due to the breakdown of reflection sym- 
metry. For r close to t„, the Hopf bifurcations con- 
tinuously split, as it is shown in Fig. |2Ia) for r = 0.1. 
This occurs in a similar way as in previous works study- 
ing the breakdown of reflection symmetry without time 
delay |l5j |. Generally, such symmetry breaking implies 
asymmetric non-simultaneous nucleations, as well as a 
complete modification of the bistability types and tran- 
sitions to synchronization. 

Figure 21 shows the region of incoherence in the t — K 
plane for two values of Alu. Due to the splitting of the de- 
generate Hopf bifurcations as r is shifted from t„, local 



extrema of incoherence appear exactly at r = r„. Ac- 
tually, the overall picture is more complex for arbitrary 
values of 7, since some peaks are washed out, and dif- 
ferent bimodal distributions may show specific shapes. 
However, in the 7 — > limit, we provide a general 
distribution-independent statement (suggested from nu- 
merical solutions of Eq. H14(l and analytical arguments): 
peaks will appear at all r = r„ values, unless there is a 
resonance (Aw/u) — p/q) what excludes n = 2mq/p, m 
even (odd) for positive (negative) K values. 

Conclusions.- We have studied the Kuramoto model 
with time delay by analyzing the linear stability prop- 
erties of the synchronized and incoherent solutions for 
bi-delta distributions Our results have been com- 

pared to previous studies on phase oscillators coupled 
via time-delayed interactions. In particular, in contrast 
to the N = 2 case |17j, the synchronous states only ad- 
mit small phase differences. Also, compared to unimodal 
distributions stability boundaries of the incoherent 
state have a completely different structure in parameter 
space, due to the presence of a new time scale in the 
system which results in a quasiperiodic pattern. Finally, 
we have considered the effect of diversity on the incoher- 
ent state exploiting the refiection symmetry of the Ku- 
ramoto model model for certain values of the time delay. 
For such values, we have achieved some analytical results 
that extend the non-delayed case studied in , and that 
help to build up a general picture for the shape of the re- 
gions of incoherence. These results may be of interest for 
the study of dynamical networks with bimodal frequency 
distribution, where the transmission times may play an 
important role on the synchronization dynamics, e.g. in 
neuroscience. For example, in the visual cortex, different 
populations of neurons form spatially separated cortical 
columns that interact with a significant time delay. In 
this context, a successful phase-model for visual process- 
ing was proposed in , but the time delay effects were 
only studied for two oscillators. 

We thank Claudio Tessone, Bernd Blasius and Jiirgen 
Kurths for fruitful discussions. E.M. was partially sup- 
ported by the European research project EmCAP (FP6- 
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